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In any valid Monte Carlo sampling that realizes microcanonical prop- 
erty we can collect statistics for a transition matrix in energy. This matrix 
is used to determine the density of states, from which most of the ther- 
modynamical averages can be calculated, including free energy and en- 
tropy. We discuss single-spin-flip algorithms, such as flat-histogram and 
equal-hit algorithms, that can be used for simulations. The flat-histogram 
algorithm realizes multicanonical ensemble. We demonstrate the use of 
the method with applications to Ising model and show its efficiency of 
search for spin-glass ground states. 

1 Introduction 

In traditional Monte Carlo sampling method, the computation of a thermody- 
namic quantity {Q) is usuaUy through a simple arithmetic average: 



where the configurations ai are generated according to a specified distribution, 
such as the Boltzmann distribution. However, it is possible to collect other in- 
formation in the same simulation, from which we can obtain better statistics, or 
estimates of quantities other than that at simulation parameters. The histogram 
method PP and multi-histogram method T collect energy histogram at a given 
temperature, from which the quantity at nearby temperature can be inferred. 
The key observation here is that the histogram of energy is related to density of 
states through H{E) oc n{E) exp{—E/kT) (in canonical ensemble). From the 
histogram, we can determine the density of states n{E). Once the density of 
states is known, we can compute most of the thermodynamic quantities at any 
temperature. 
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Histogram method has been found to be an excehent tool for study criti- 
cal phenomena. Further improvement can be made by collecting 'high-order' 
statistics, i.e., the transition matrix With histogram method, each 

configuration provides just an '1' to a histogram entry, while in transition ma- 
trix method, each configuration gives several numbers of magnitude about A'^ 
to the transition matrix elements, thus variance reduction is expected. One 
of the most appealing feature of transition matrix Monte Carlo is an easy and 
straightforward way to combine several simulations. Additionally, we can use 
any valid sampling algorithm in a generalized ensemble which realizes the mi- 
crocanonical property that states with same energy have the same probability. 
The flat-histogram algorithm [3 |S] is such an algorithm that realizes multi- 
canonical ensemble in which the energy histogram distribution is a constant. In 
the following, we present the transition matrix Monte Carlo method, introduce 
the flat-histogram and other related algorithms. We discuss the performance 
of algorithms with examples from Monte Carlo simulation results of the Ising 
models. We summarize in the last section. 

2 Transition Matrix Monte Carlo Method 

First, we give the definition for the transition matrix. Let W{(j a') be the 
transition probability of the states from a to a' of a Markov chain. To be 
definite, we consider a single-spin-flip dynamics with a canonical distribution, 
but the formalism is general. The transition matrix in the space of energy from 
E to E' is 

^ ^ m^x-a'), (2) 

^ ' E(<y)=E E{(j')=E' 

where the summations are over all initial states a with energy E and all final 
states fj' with energy E' . Estimates of the transition matrix can be obtained 
during a Monte Carlo sampling, where the summation over E divided by n{E) 
is interpreted as a microcanonical average of the state-space transition proba- 
bilities of the Markov chain, i.e., 

T{E^E')= J2 {W{a^a'))^, (3) 

E{a')=E' 

The expression can be further simplified if we consider single-spin-flip dynam- 
ics, with a spin choosing at random. In this case, '^E{a')=E' ^{'^ ~* '^') — 
j^N{a, AE)a{E E'), where N{a,AE) is the number of sites such that a 
spin-flip causes the energy increasing by AE = E' — E in the current state 
a. It is also the number of possible moves that one can make to change the 
energy by AE. Note that J2ae -^('^i = -^j where TV is the number of 
sites. A common choice of the single-spin-flip rate a{E E') is the Metropolis 
rate min(l, exp(— Ai?/fcT)) . Since this factor is a function of E and E' , the 
microcanonical average {■■■)£ is performed over N(a, AE) only. We have 

T{E ^ E') = ^{N{a, AE))^a{E ^ E') = T^{E ^ E')a{E ^ E'), (4) 
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where we have defined a normalized N{a, AE) as the infinite temperature tran- 
sition matrix. 

The eigenvector corresponding to the eigenvalue 1 of the transition matrix 
is the probability of finding states with energy E. It is also proportional to the 
histogram H{E). To determine the density of states, a numerically better choice 
is from the detailed balance. This gives us the relationship between histogram 
and transition matrix: 

H{E)T{E E') = H{E')T{E' -> E). (5) 

If we use the fact that H{E) cx n{E) exp{—E/kT) and Eq. |@J, we obtain the 
so-called broad- histogram equation [U] ^1 ^2 ^] 

niE){Nia, E' - E)) ^ = n{E'){N{a', E - E'))^,. (6) 

This is one of the basic equation for determining the density of states, as well 
as for the flat-histogram algorithm below. 



3 Flat-Histogram Algorithm 

Any sampling algorithm that can realize microcanonical property, i.e., the distri- 
bution of the states is a function of energy only, can be used to collect statistics 
for {N{a,AE))E- Using a canonical ensemble simulation, we need dozen tem- 
peratures in order to cover all the relevant energies. However, comparing to 
multi-histogram methods, the combination of data at different temperatures is 
very easy, we simply add up the matrix elements and then properly normalize. 

Multicanonical ensemble ^H] is a particularly good choice for the collection of 
transition matrix elements, since it reaches all energy levels with equal probabil- 
ity. Multicanonical ensemble is defined to be H{E) — const, or the probability 
of configuration P(cr) oc l/n(^E{a)) . It is purely an artificial ensemble designed 
for computational efficiency. To realize the multicanonical ensemble, we can per- 
form a single-spin flip with a flip rate of min {l,n{E)/n{E')). However, smce 
the density of states n{E) is not known beforehand, we have proposed to use 
the count number N{a, AE). From the broad-histogram equation, Eq. the 
ratio of 7i(-) is related to the ratio of N(-), we have 

f '\ ■ (^ {N{a',E~E'))E' \ 
°^"^"^ = "^^H ' {N{a,E'-E))J - 

This is our flat-histogram flip rate. Although the microcanonical average {N{a, AE)) e 
is also not available before the simulation, it can be obtained approximately dur- 
ing a simulation. We use the instantaneous value and running average to replace 
the exact microcanonical average. Numerical tests have shown that this proce- 
dure converges to the correct ensemble for sufficiently long runs. For realizing 
truly a Markov chain, it is sufficient for a two-pass algorithm. The first pass is 
as before. In the second pass, we use a multicanonical sampling rate, using the 
density of states determined from the first pass. 
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A variation of the algorithm is equal-hit algorithm which combines the N- 
fold way method with a flip rate that gives an extended ensemble that is 
uniform in probability of visiting each new energy. Reference [5] gives more 
extensive discussions, as well as comparison with Wang-Landau method |15| . 

4 Some Results 

As there are more detailed balance equations among the transitions of different 
energies than the number of energy levels, we determine the density of states 
from the transition matrix by solving a least-squares problem, or more generally, 
a nonlinear optimization problem. The optimization can be done either in the 
density of states n{E), or in the transition matrix elements T{E E'). There 
are a number of constraints that the transition matrix must satisfy. The trivial 
one is the normalization, J^AE^^i^ E + AE) = 1. There exists a rather 
interesting constraint, known as TTT identity, as well: 

T{E -> E')T{E' E")T{E" -> E) = T{E -> E")T{E" E')T{E' E). 

(8) 

These constraints complicate the optimization algorithms. 

While any of those extended ensemble methods reduce their efficiency as the 
system size increases, the accuracy of a two-pass flat-histogram/multicanonical 
simulation is rather good for a given fixed amount of CPU times The 
method can also give excellent result for large systems, such as a 256 x 256 Ising 
lattice jJT], using a parallelized version of the program. The method is also 
applied to a lattice protein model, the HP model, with good performance 

A possible measure of computational efficiency is through the tunneling 
times. The tunneling time is defined to be the number of Monte Carlo steps 
in units of a lattice sweep {N basic moves), for system making a pass from the 
highest energy level to lowest level, or vice versa. For the two-dimensional Ising 
model, this tunneling time diverges with system sizes according to L^'^ which 
is worse than standard random walk. On the other hand, for spin-glasses with 
complicated low-temperature free-energy landscape, the tunneling time is much 
larger. It is about L'*'^ (19) in two dimensions and L^'® ^5 in three dimensions. 
Another measure for spin glasses is given by the average first-passage times. It 
is defined as the average number of sweeps needed to reach a ground state. It is 
found in ref. (20| that the first-passage time diverges exponentially rather than 
according to a power. In any case, the equal-hit algorithm performs comparable 
to 'extremal optimization' '21' which is an optimization algorithm inspired from 
self-organized criticality. 

5 Conclusion 

By collection the transition matrix, more information is obtained about the sys- 
tem, giving more accurate results. The effect of using transition matrix is more 
dramatic for small systems. Although the transition matrix analysis of data can 
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be used with any simulation algorithms, extcnded-cnscmble-bascd algorithms, 
such as flat histogram algorithm, are excellent choices. The eflSciency of the 
flat-histogram related algorithms has be studied. 
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